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Non-normal transient growth of disturbances is considered an essential prerequisite for 
subcritical transition in shear flows, i.e. transition to turbulence despite linear stability 
of the laminar flow. In this work we present numerical and analytical computations of 
transient growth covering all linearly stable regimes of Taylor-Couette flow. Our numer- 
ical experiments reveal comparable energy amplifications in the different regimes. For 
high shear Reynolds numbers Re the optimal transient energy growth always follows a 
Re 3 -scaling, which allows for large amplifications even in regimes where the presence of 
turbulence remains debated. In co-rotating Rayleigh-stable flows the optimal perturba- 
tions become increasingly two-dimensional as the optimal axial wavenumber goes to zero. 
In this limit of axially invariant perturbations we show that linear stability and tran- 
Mh sient growth are independent of the cylinders' rotation-ratio and we derive a universal 

i?e 3 -scaling of optimal energy growth using WKB-theory. Based on this a semi-empirical 

t-H formula for the estimation of transient growth valid in all regimes is obtained. 
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Q 1. Introduction 

The flow of viscous fluid between two coaxial independently and uniformly rotating 
cylinders, Taylor-Couette flow, is a paradigmatic system to study the stability and dy- 
namics of rotating shear flows. For simplicity we assume here that the system is infinite 
in the axial direction so that the annular geometry is uniquely determined by the dimen- 
$Jj sionless radii ratio r\ of the inner and outer cylinder. 

The corresponding laminar Couette flow is determined by the inner and outer Reynolds 
numbers Re^ and Re which are proportional to the rotation frequencies of the respective 
cylinders f2i and Q . It is well-known that the base flow's stability does not only depend 
on the shear magnitudes of Re-i and Re but also changes qualitatively with their ratio. In 
particular, inviscid Couette flow is linearly stable if and only if the fluid particles' angular 
momentum increases in radial direction according to Rayleigh's criterion ( |Rayleigh|1917[ ). 
Consequently, inviscid stability solely depends on the ratio Rei/Re in this case. For 
viscous fluids this leads to a complex interplay of shear- and centrifugal stability which 
governs the competition between laminarity and turbulence. 
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Figure 1. T aylor-Couette flow regimes in the Rei-Re plane by the parametrization of |Dubrulle| 
et aZ.| ( 2005 1 (n = 0.5). The rotation number Rq uniquely determines the regimes I to IV, whereas 
the shear Reynolds number gives the magnitudes of Rei and Re as visualized in the plot 
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Table 1 . Parametrization of the Taylor-Couette flow regimes by the rotation number Rn as 
visualized in figure [l] lines of constant Rq are axes meeting in the origin 



In an attempt to separate these effects we adopt the parametrization by |Dubrulle| 
et al. ( J2005 1 using shear Reynolds number Re and rotation number Rq to parametrize 
the Rei-Re plane (see figure [l]). As the name suggests Re ~ fl — fli is a measure for 
the (absolute) shear in the flow, whereas Rq depends solely on the ratio fl /H,i. 

For viscous fluids Rayleigh's criterion then yields certain ranges of Rq in which the 
base flow remains linearly stable up to arbitrary large Re. These Rayleigh-stable regimes 
correspond to I and II in figure [l] The remaining regimes III and IV are Ray leigh-unst able 
so that the laminar flow develops linear instabilities for Re — > oo. In practice these appear 
already at moderate Re — O(10 3 ) except when approaching the boundaries of regimes I 
and II ( Taylor|[i"923 ) . Table fT] gives the parametrization of the various regimes via Rq . 

We subdivide the Rayleigh-stable regime according to the base flow's angular velocity 
profile B : The quasi- Keplerian regime II is characterized by d r Q B < 0, i.e. radially 
decreasing angular velocity, whereas regime I is defined by a positive gradient d r ft B > 0. 
In figure [T] these domains are separated by the solid body line given by Rei = r\Re . 
Due to the absence of shear for these configurations fi B is constant, corresponding to 
a flow profile that equals a rotating solid body. The transition from regime II to III 
defines the Rayleigh line where Rayleigh's stability criterion ceases to be fulfilled and 
centrifugal (linear) instability sets in. In experiments this results in the formation of a 
new stationary flow, characterized by the famous toroidal Taylor-vortices ( Taylor|1923 l. 
Similar instabilities and associated patterns occur in a wide range of the counter-rotating 
regime IV which is characterized by opposite signs of Rei and Re . For these flows very 
good agreement between linear stability analysis and experimentally observed stability 
boundaries has been achieved for moderate Re. 

However, similar to plane Couette and Poiseuille flow (compare Romanov|1973 Davey 



1973 Drazin & Reid 1981 1 certain Taylor-Couette flows may undergo subcritical tran- 



sition to turbulence in absence of unstable eigenvalues. This phenomenon has been ob- 
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served both in the Rayleigh- unstable regime IV ( |Coles|1965 l as well as by | Wendt (19331 
and Taylor ( 1936[ ) at the lower boundary of the Rayleigh-stable regime I (i.e. for fixed 
inner cylinder see figure [I]) . Recent studies by Borrero-Echeverry fc Schatzj ( J2010 ) have 
confirmed the rapid lifetime increase of turbulent spots with the Reynolds number in the 
latter setting. Hence, we may infer the existence of subcritical turbulence within regime 
I in spite of the lack of experimental and numerical data for such flows. 

On the other hand, the existence of turbulence remains controversial in the equally 
Rayleigh-stable quasi-Keplerian regime II ( Yecko|20"04 Ji et «/.|2006 Paoletti & Lathrop 



2011 Balbus 2011). As the name suggests these flows are of great importance in mod- 



elling astrophysical objects with Keplerian velocity profiles such as accretion disks (for 
details see Pringle 1981). However, endcap effects render this regime difficult to explore 
experimentally. In fact, Avila (2012) has shown state-of-art Taylor-Couette apparatuses 
unsuited to adequately produce the respective flow fields at the required Reynolds num- 
bers. Based on i?e-bounds derived from a variational formulation of the stability problem, 



Busse (2007) conjectured that turbulence cannot exist in the quasi-Keplerian regime. Yet, 



this result is predicated on the hypothesis that the extremalizing vector fields are inde- 
pendent of the streamwise coordinate. To the best of our knowledge there is no general 
proof ruling out the existence of turbulence in the literature. 

Regardless whether linear or nonlinear, stability analysis boils down to the evolution 
of initial perturbations to a stationary state. For stationary flows, the development of 
the perturbations' energy is given by the Reynolds-Orr equation which is valid both 
for fully- nonlinear and linearized dynamics (Schmid & Henningson 2001). Remarkably, 
this implies that nonlinear instabilities may exist only if the linearized Navier-Stokes 
equations have solutions that grow in energy, i.e. transition requires linear growth. 

At first glance, this theory seems contradictory to subcritical transition being a mani- 
festation of nonlinear instability despite linear stability. However, the paradox is resolved 
by the non-normality of the linearized Navier-Stokes operator, i.e. the non-orthogonality 
of its eigenmodes (Kato 1995). This potentially allows for transient growth of infinitesi- 
mal perturbations (Boberg & Brosa 1988; Trefethen et al. 1993), i.e. temporary energy 
growth even in case of linear stability (as illustrated e.g. by Grossmann (2000)). Like 
in other flow geometries the non-normality of the Taylor-Couette operator grows with 
the shear Reynolds number Re so that the maximum energy amplification, G max , may 
reach several orders of magnitude at sufficiently large Re (Reddy & Henningson 1993). 
For instance numerical simulations by Yecko ( 2004 ) of the rotating plane Couette geom- 
etry showed an asymptotic scaling of G max ~ Re's for one configuration of the highly 
controversial quasi-Keplerian flows in the limit Re — > oo. 

For counter-rotating Taylor-Couette flows Meseguer ( |2002 ) has studied optimal tran- 
sient growth at the subcritical stability boundary Rex{Rn) measured by Coles (1965). 
Most prominently, he partly observes a strong correlation and finds a sharp threshold 
value Gmax.T = 71.58 ± 0.16 for relaminarization in the experiments. These results rein- 
force the potential significance of non-normal growth in subcritical transition. 

This article is concerned with transient growth in all regimes of linearly stable Taylor- 
Couette flows identifying universal properties, especially in the limit of high Reynolds 
numbers. After briefly presenting the governing equations of the Taylor-Couette problem 
and our numerical formulation in <|2] and <|3] we discuss some tests of the method and 
numerical issues of transient growth computations in <HJ In SJ5]the main numerical results 
for the asymptotic scaling G max ~ Re a of optimal transient growth and the correspond- 
ing optimal perturbations are presented. Furthermore, a semi-empirical formula for the 
estimation of G max by Re and the cylinder radii ratio rj is obtained. The latter is revealed 



4 S. Maretzke, B. Hof and M. Avila 

to be universal by the analytical results for axially independent perturbations derived in 
g6J For such disturbances we further verify the characteristic scaling G max ~ i?e s via a 
WKB-approximation to the linearized evolution equations in <|7] In the final section SJ8] 
we discuss our results and draw some conclusions concerning subcritical instability. 



2. The linearized Taylor Couette problem 

2.1. Principal equations 

We consider an incompressible Newtonian fluid with kinematic viscosity v confined be- 
tween two coaxial independently rotating cylinders with radii r[ < r' a that are infinite 
in the axial direction. Nondimensionalized with the gap width d := r' a — r[ as length 
scale, viscous time v~ x d 2 and the pressure scale v~ 2 d 2 the system is governed by the 
dimensionless incompressible Navier-Stokes-equations and continuity equation 

d t v = -(v V)v- Vp + Av (2.1a) 

V-v = (2.1b) 

where p is the reduced pressure and v the velocity field of the fluid. 

The independent variables are the viscous time t and the spatial vector x parametrized 
in cylindric coordinates x =: (r, ip, z) 1 . The dimensionless geometry parameters are given 
by ri := r\d~ , r := r' d~ l and the radii ratio r\ := r^r" 1 . Let fli and fl be the (constant) 
angular velocities of the inner resp. outer cylinder. Defining the inner and outer Reynolds 
numbers Ret := -r^fli and Re :— -r' il the no-slip boundary conditions at the inner 
and outer cylinder walls read 

V\r=rt = R e i e ip an d t, |r=r = Re £<p (2-2) 

where e r =: (1, 0, 0) T , e v =: (0, 1, 0) T and e z =: (0, 0, 1) T denote the orthonormal radial, 
azimuthal and axial unit vectors. 



A well-known solution of the boundary value problem ( 2.1 ) and (2.2 1 is laminar Couette 
flow (v B ,p B ), given by 



B\ 1 B 2 

v B =[Ar+ — )e v and p B = -A 2 r 2 + 2ABln(r) (2.3a) 

r ) 2 2r 2 

l + n (1-77)(1 -77 2 ) V ! 



In order to investigate its stability the equations (2.1) are linearized about (v B ,p B ) 
yielding the linearized Navier-Stokes-equations for the evolution of an infinitesimally 
small perturbation (u,q): 

d t u = -(v B ■ V)u - {u ■ V)v B -Vq + Au (2.4a) 

V-m = (2.46) 

ui r=ri = u\ r=ro = (2-4c) 

By a Fourier ansatz in the azimuthal and axial coordinates u(r, <p, z) := u(r)e 1 ( nip+kz \ 
q(r, ip, z) := q( r )eA nip+kz,> for k £ K, n S Z the evolution equation can be written as 

d t u = Cu- V r q. (2.5) 

Herein a subscript r for an operator T denotes the conjugate with e l ( n( f+ kz ) ; i. e . % := 
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e -i{ nv +kz) 1 - ei (n V +kz)_ The operator £ is given by (JMeseguer| |2002| ) 



Cu = -{y B ■ V) r u - (u ■ V)v B + A r u 




(2.6a) 



c — c — c 



= v+v- 



n 2 + l 



-k 2 - 



2 R 2in 
^ r v> r > r 2 



r 
d r and 2? + 



2in 



V+vt 



(2.6b) 



d r + -. The domain of admissible velocity 



with the abbreviations V 

fields u = (u r , u v , u z ) J in equation (2.5) is the twice continuously differentiable subspacc 

V:={t)£i 3 n ^ 2 ((r 2 ; r )) : v{n) = v(r ) = 0, V r • v = 0} (2.7) 

of the Hilbert space H 3 . Here we define H := h 2 ((ri; r )) with the inner product 



(v) 



H • ' 



-; (91,92) i-» / qlq 2 rdr. 



(2.8) 



where the superscript * denotes the conjugate transpose of a scalar, vector or matrix. The 
square of the induced norm ||m|| 2 of (•, •) := (•, -) H3 is proportional to the total kinetic 
energy of the perturbation u and is thus called the energy norm. 

A modal ansatz in the time coordinate t, i.e. u := 
yields the eigenvalue problem 



U\e and q := q\e M for A G 



Att A = Cu x - V r q\, (u x ,q\) E V x 



(2.9) 



For the axisymmetric case n = DiPrima & Habetler ( 1969 ) have shown the discreteness 



of the eigenvalues {A} and completeness of the corresponding generalized eigenfunctions 
in V. If we assume that this remains true for n 7^ then the laminar Couette flow (2.3) 



is linearly stable if and only if all eigenvalues of (2.9) have negative real parts 



2.2. The parameter space for transient growth 



In addition to the geometric parameters Rei, Re and n the evolution problem (2.5) 



depends on the azimuthal and axial wavenumbers n and k. Due to the cylindric symmetry 
of the Taylor-Couette geometry the parametric analysis may be confined to Rei, n,k ^ 



(see Meseguer & Marques (2000) for details). The parameter 77 e (0; 1) determines the 
curvature of the system and thus the rotational influence. The limit -q — > 1 corresponds 



to plane Couette flow as demonstrated by Hristova et al. ( 2002 ) with respect to transient 
growth, whereas 77 — >- implies infinite curvature at the inner cylinder wall. 

For reasons discussed in §1 we introduce the shear Reynolds number Re and rotation 
number Rq. Assuming Rei ^ and Rei ^ r]Re the mapping (Rei,Re ) <-> (Re,R&) is 
one-to-one so that the flow parameters A, B can be expressed via Re and Rq: 



Re 



2\nRe -Rej\ 
1 + ry 



and Rq 



(l-ri)(Rei + Re ) 



rjRe — Re^ 



sgn(R n )Re 
A = 7, ( H 



1) and B = - 



sgn(Rn)yRe 
2(1-77)2 ■ 



(2.10a) 



(2.106) 



S. Maretzke, B. Hof and M. Avila 



Computing the commutator of the operator C given by (2.6) with its adjoint C* , 
[£* , £} — 0(Re ), reveals its non-normality scaling with the shear Reynolds number. The 



eigenspaces are therefore non-orthogonal to one another ( Kato|[l995 1 , which potentially 



allows for significant transient growth at large Re. Detailed discussions of this mechanism 
can be found in ( ]Grossmann||2000[ ) and (Sch mrd fc Hennmgson|2001| pp. 99-101). 

As a consequence initial perturbations u(0) may temporarily grow in energy before 
they ultimately decay - even if C has only stable eigenvalues AeC with Re(A) < 0. The 
maximum transient growth at time t ^ is given by G(t) :— sup|U( )||=i ll M WI| 2 - If f ne 
evolution of u is linear, i.e. d t u = Au where A a linear operator G can thus be expressed 
using the operator norm (Trefethen et aZ.||l993 ): 



G{t):= sup ||m(£)|| 2 = sup || exp(„4i)it(0)|| 2 = || exp(„4f)|| 2 (2.11) 

||u(0)||=l ||u(0)||=l 

If || • || denotes the energy norm G(t) is equal to the greatest kinetic energy amplification 
an initial perturbation w(0) € V can attain. 

For a Taylor-Couette flow configuration given by the parameters Rei, Re and r\ the 
optimal transient growth is defined by G max := sup t n k G{t). A perturbation u with 
||m(0)|| = 1 is called optimal if ||ti(£)|| 2 = G max for some t ^ 0. Note that G max is finite 
if and only if all eigenvalues of C are stable. 



3. Numerical Formulation and Implementation 

3.1. The Galerkin method 



The eigenvalue problem (2.9) is numerically solved using a Galerkin method. The imple- 
mentation is similar to the Petrov-Galerkin method described by Meseguer & Marques 
(2000) and Meseguer et al. ( 2007[ ), however based on Legendre rather than Chebyshev 



polynomials so that trial and projection basis are identical. 

The basis choice is U := {w^I^n where v} m and u^ are defined according to table 
[2] for different wavenumbers n, k. The functions h m and g m are given by 

h m (r) := r(l - x 2 )L m (x) and g m (r) := r(l - x 2 ) 2 L m (x) for re[r,;rj. (3.1) 

where L m is the Legendre polynomial of degree m and x := 2r — (1 + 77) (1 — rf) . Then 
every v? m satisfies both the continuity condition V r -M^ = and the boundary conditions 
by definition since by construction 

h m (n) = h m (r ) = g m (n) = g m (r ) = s^(n) = g' m {r ) = 0. (3.2) 



The problem is discretized by truncating U at the polynomial resolution N G N, i.e. 

i expanding possil 

m<N,j = l,2 u 'm Ul m 



defining Un '■= { u mYm<N > anc I expanding possible solutions U\ to the eigenvalue problem 
(2.9) in terms of U N , u x := J2 m <N 7=1 2 a m u m- Plugging tms ansatz into equation (2.9) 



and projecting on some u\ yields 

A Y, H<> < = E < u ?' £ <) < - H V r g) . (3.3) 

m<N m<N V v ' 

3=1,2 j=l,2 =0 

The pressure terms vanish due to the boundary- and divergence conditions. 



Thus, the equations (3.3 1 for alii < N, i = 1, 2 can be written in form of a 27V x 27V 



generalized eigenvalue problem 

\Ga = Ha with G :=(«<», H :=(«£<» (3.4) 
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0, fc = 



0, fc/0 n/0, fc = 0n/0, fc/0 




' —ikrg m 
D+(rg m ) 



, /'/. 




Table 2. Spe ctra l basis function s for m £ N us ed for the dis cretization of th e eigen value 
problem (2.9 1 via equations (3.1l and (3.3 1 according to Meseguer et al. (20071 



for the coefficient vector a := (aj, • 

matrices G being Hermitian positive definite (Meseguer & Marques 2000). 



*jV-l!"0' • • 



i Oat_i) t where G and H are 2N x 27V- 



3.2. Computation of transient growth 

Now let Q := {q 1 ,...q 2 Ar} be the eigenf unctions corresponding to the eigenvalues 
A := {Ai,...A2at} and eigen(-coefficient-)vectors {ai, . . . cl2n} solving the generalized 



^2N 



eigenvalue problem (3.4 1. Consider some perturbation expanded in Q, i.e. u — J2i=i biQi 
where b — (bi, . . . , &2Af) T denotes the time-dependent coefficient vector. Since the q i are 
(approximate) solutions to the eigenvalue problem (2.9) it follows that 



b(t) = exp (diag(A)t) 6(0) 



(3.5) 



where diag(A) denotes the diagonal matrix constructed from A and exp is the matrix 
exponential. Thus the evolution of the perturbations kinetic energy reads 



2 A 



Ml 2 = (u,u) = Y, b*b j (q i ,q j ) = b*Mb=\\Fbf 2 = ||Fcxp(diag(A)i) b(0)\\$ (3.6) 



Here M is the Hermitian positive definite Gramian matrix M := ((gj,q,-)), M = F* F 
a Cholesky decomposition and || • |J2 denotes the standard 2-norm on C 2N . Hence, the 
maximum transient growth at time t ^ is given by (see |Meseguer||2002 1 



G{t) = sup ||it(£)|| 2 = sup ||Fexp(diag(A)i)b|| 



*(o)|| 



\Fb\\ 



-Fb 



sup ||Fexp (diag(A)t) F l v\\\ 

|w||2 = l 



| F exp (diag(A)t) F 



-1||2 

lb- 



(3.7) 



So G(t) is equal to the squared maximum singular value <r 2 of Fexp (diag(At)) F . 
Moreover, if Vq denotes the corresponding right-singular vector Fvq is the initial Q- 
coefficient vector of a perturbation which attains optimal transient growth at time t. By 
means of singular value decomposition we thus obtain both maximum transient growth 
G(t) and corresponding perturbations in the finite-dimensional subspace spanned by Q. 
This yields a lower bound to the maximum attainable by arbitrary initial conditions in 
V. As discussed in §4.2|we find convergence of this estimate to the total maximum. 
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3.3. Outline of the code 

By definition of Un only integrals over polynomial functions have to be evaluated in 
order to calculate G and H. Hence, these are computed exactly using Gauss-Legendre- 
quadrature with Gauss-Lobatto collocation points of degree M where M > N + 6 (see 
69 ff.). 



Canuto et al. 2006 



pp. t>y it.). Moreover, the derivatives in the operator C are imple- 
mented by means of the corresponding differentiation matrices given in |Canuto et al.\ 
p006| p. 76). 



The Legendre polynomials are determined using recursion formulas starting at the 
corresponding Legendre-Gauss-Lobatto collocation points. The latter are computed via 
Newton's method starting from the respective C/ie&j/s/iew-Gauss-Lobatto points as initial 

pp. 69 ff.). 



guess (Canuto et al. 2006 pp. 69 ff.). The iteration is performed up to a maximum 
deviation of e = 10 _ib between two subsequent iterates. 

The code used in this work is based on the scientific computing package Scipy for the 
interactive language Python. The linear algebra algorithms are provided by the package 
Scipy. Linalg based on the standard ATLAS, LAPACK and BLAS implementations. 

The optimization of G in the time coordinate £ € [0; £ C ut] and in the continuous wave 
number k £ [0; fc cu t] are performed via the Scipy .Optimize implementation of Brent's 
method (for details see Press et aT||2007 section 9.3). With respect to the discrete wave 
number n £ {0, 1, ... , n C ut} G is optimized by brute force. If the optimal transient growth 
is found at the upper boundary of the considered domains, i.e. for t — £ cut , k — k cnt or 
n = "cut the respective intervals are enlarged in subsequent steps until a local maximum 
is located in their interior. 



4. Numerical issues 

In this section the performance of the numerical implementation presented in section 
[3] is tested by comparison to results in literature. Furthermore, eigenvalue and transient 
growth convergence are studied for test cases in order to justify the choice of polyno- 
mial resolution N used to obtain the numerical results in f|5] We find that the optimal 
transient growth may converge without the Y-shaped spectrum being properly resolved. 
This observation suggests a lesser significance of the spectrum in transient growth com- 



putations contravening elaborations by IReddy & Henningson ( 1993 ) for channel flows 



4.1. Eigenvalue Decomposition 



Our discretization of the eigenvalue problem (2.9) has been tested against the results 



on eigenvalue-critical Reynolds numbers presented in Krueger et al. ( 1966 table 2) as 
well as by replication of the plotted spectra given by Gebhardt & Grossmann ( 1993| fig. 
3a-d) . Agreement within the respective accuracies has been found. Additionally, we have 
compared our Galerkin method to the Petrov-Galcrkin scheme of Meseguer et al. (2007). 



No significant deviations are found between the converged spectra. 

For these methods we study the convergence of the approximated least stable eigen- 
value Aj^ as the number of Legendre - resp. Chebyshev polynomials N is increased. In 
figurebkhe relative errors |Af — A x '"'HAj rcf | _1 compared to (converged) reference values 
A-l rof are plotted against N. The test parameters are Rq = —2, r\ = 0.5, n — 5 and k = 1 



at shear Re ynold s nu mber s Re = 8000 (figure |2 (a) [ > resp. Re = 128000 ( |2(b)[ ). 

The plots |2(a) and |2(b) show plateaus of non-convergence for low N which are due to 
the difficulty of identifying the resp. eigenvalue in a non-converged spectrum. For moder- 
ate N € [20; 45] (fig. |2(a)[ ) resp. N £ [45; 90] (fig. |2(b)[ ) spectral accuracy, i.e. exponential 
convergence rates, is observed for both methods. Notably however, the convergence turns 



Transient growth in Taylor- Couette flows 



Galerkin 
?e:-ov-Ga e'< ■ 



Polynomial resolution TV 

(a) Re = 8000 



Galerkin 
3 e:'ovGa t--< ■ 



Polynomial resolution N 

(b) Re = 128000 



Figure 2. Convergence of the least stable eigenvalue Ai of C for Rn — —2, r\ — 0.5, n — 5, k — 1 
and ije = 8000 (a)), Re — 128000 (b)) computed u sing our Galerkin method (triangles) and the 
Petrov- Galerkin scheme of Meseguer et al. (20071 (crosses); Ai denotes the approximation to 



Ai computed using TV Legendre resp. Chebyshev polynomials; |Ai — A x ref |/|Aj rct | is the relative 
deviation of A^ from the converged result X x rBf 



Meseguer (2002) Present work (TV = 50) 



tlG-i ^t6o "'max %ax *- r : 



max ' 'max a max 



"'max ^~ r n 



591 -2588 

523 -2975 

473 -3213 

405 -3510 



10 

11 
11 
11 



1.994 
1.996 
1.920 
1.839 



71.36 
71.58 
71.64 
71.75 



10 

11 
11 
11 



1.997 
1.998 
1.922 
1.841 



71.58 
71.81 
71.87 
71.99 



Table 3. Optimal transient growth G max := sup nkt G(t) according to Meseguer (2002 
1) and own results; parameters are r\ — 0.881 and TV — 50; n 
resp. axial wavenumbers which attain optimal transient growth G 



Table 



out to be significantly quicker in case of the Legendre polynomial based Galerkin method 
presented in this work: spectral accuracy is attained using significantly fewer polynomi- 
als and the limiting machine precision is reached already for TV = 43 (Re = 8000) and 
TV = 83 (Re = 128000) compared to TV = 62 resp. TV = 104 in case of Petrov-Galerkin 
scheme (see figure [2]). 

The required resolution TV for convergence grows with the shear Reynolds number 
Re and - much more significantly - as soon as subsequent, more stable eigenvalues are 
considered. In fact, it turns out to be numerically impossible to resolve significant parts of 
the eigenvalue spectrum for Re ^ O(10 5 ). This also affects the computation of transient 
growth discussed in the next subsection. 



4.2. Computation of Transient Growth 

In tableppur results on the optimal transient growth G max := sup„ k t G(t) for n = 0.881 
and the corresponding optimal wave numbers n max are compared to the numerical data 
of Meseguer] ( 2002 Table 1). Beyond a deviation ^ 0.3% in fc max and G max the results 



are in perfect agreement to one another. 

The convergence of the maximum transient growth G shows remarkable characteris- 
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Figure 3. Eigenvalues \i (top), time dependent maximum transient growth G(t) (middle) and 
modulus |w r (r)|, I^MI an d l w z( r )l of the radial, azimuthal and axial components of the pertur- 
bation u attaining optimal growth sup t>0 G(t) (bottom) approximated by different resolutions 
TV; parameters: Re = 10000, Rn — — 2, r\ = 0.8, n = 5 and k — 1 



tics which partly contravene the significance of the linearized Navier-Stokes operator's 



spectrum for such computations claimed e.g. by Reddy & Henningson ( 1993 ) 



These features are discussed with reference to the example displayed in figure [3} for 
three different resolutions N £ {5, 15, 50} (corresponding to figures 3(a) 3(b) and 3(c) ) 
the eigenvalues (top), the evolution of the maximum transient growth G(t) (middle) and 
the components' moduli |u r |, \u v \ and \u z \ of the corresponding optimal perturbation 
tx(0) are plotted for comparison. The example parameters are Re = 10000, Rq = —2.0, 
i] = 0.8, n = 5 and k = 1. 

A few aspects are noteworthy: around its maximum G is already surprisingly well 
approximated by only N = 5 Legendre polynomials whereas the optimal perturbation is 
far from its actual shape (see figure [3(a)] ) . For N = 15 (figure [3~(b)| the curve {(£, G(t))} 
is converged within an error sj 1% while its maximum is even approximated up to « 
0.01%. Likewise, the optimal perturbation ii(0) is practically converged. At the same 



time the characteristic Y-like structure of the eigenvalue spectrum (compare Gebhardt 



|&: Grossmann 1993) is by no means well resolved for N — 15 not to mention N = 5 
(top). In fact, it takes as many as N = 50 polynomials for convergence of the two 



meeting branches (see figure 3(c)). However, this does not seem to affect the transient 
growth quantities - even though the converged spectrum in figure 3(c)| (top) is even much 
more stable as a whole than its approximation for N — 15 (figure ref fig 5.2.1b). 



In contrast to these observations Reddy & Henningson ( 1993 1 stress the significance 



of the two eigenvalue branches and especially their meeting point for transient growth 
in channel flows. As for Taylor-Couette flow this is only confirmed if the Y-structure is 
resolved in the first place. This turns out not to be necessary which is a lucky circumstance 
in two respects: on the one hand the two branches consist of 0(Re a ) discrete eigenvalues 
for a « I rendering their convergence numerically infeasible for Re ^ O(10 5 ). This 
reflects the operator's continuous spectrum in the limit Re — > oo. 



On the other hand the standard Cholesky decomposition of the matrix M (see [3.2 1 
tends to fail at large Re if the eigenvalue spectrum is over-resolved. In the example 
shown in figure [3] this happens for N > 51 - just as the crucial meeting point is resolved. 
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Maximum Re 8000 16000 32000 64000 128000 256000 512000 1024000 2048000 

Resolution N Re 31 38 47 58 71 88 107 131 159 

Table 4. Canonical resolutions Nr £ for the computation of optimal transient growth G max for 
Re below the given upper bounds; determined by the convergence of G max for r\ = 0.2 



Accordingly, one might expect to miss a sudden jump in the maximum transient growth G 
if the method breaks down precisely at this point. Note however that no such is observed 
in those cases where the meeting point may still be resolved, i.e. for smaller Re. 

We may thus conclude that the transient growth of the linearized Taylor-Couette 
operator C is already converged while its approximated spectrum is still far from its 
natural shape. Startling at first glance, this is yet another manifestation of transient 
growth's non-modal nature: the non-eigendirections are those of significance. 

Nevertheless, numerical artifacts in form of spurious unstable eigenvalues have to be 
avoided by choosing sufficiently high resolutions N . However, N must not be too large 
either in order to keep the Cholesky decomposition stable (although preconditioning or 



more stable algorithms such as the one presented by Ogita & Oishi (2012) might be 
another alternative). For a given set of parameters T), Re, Rq it turns out that resolving 
the transient growth peak for optimal wavenumbers n = w max , k = fc max tends to require 
the highest resolutions. Moreover, the necessary N are widely independent of Rq and at 
least of the same magnitude for different n. Here greater curvature, i.e. n — > 0, results in 
slower convergence. Consequently, for practical computations suitable resolutions N Re are 
determined for different ranges of Re by the convergence of (computationally challenging) 
test cases, more precisely less than 0.3% deviation in the optimal transient growth for 
7? = 0.2andiVe [N Re - 3; N Re }. 

For greater n lower resolutions N may be sufficient and greater Reynolds numbers than 
Re = 2048000 might be resolvable. However, universal convergence for any parameters 
i?o, r), n and k within about 1% may be assumed if N is chosen according to the resulting 
canonical resolutions N Re given in table HI They are found to approximately follow a 
power law of the form N Re = N Re a with N = 2.28 ± 0.06 and a = 0.293 ± 0.002. 

Starting from these, N is temporarily reduced in subsequent steps whenever the Cho- 
lesky decomposition fails and temporarily increased if unstable eigenvalues occur in order 
to identify possible numerical artifacts. In case of converged unstable eigenvalues the 
computation of the matrix M and thus of the transient growth is confined to the stable 
eigenmodes in Q in agreement with the analysis oflMeseguer (2002). 



These computation guidelines have been applied to obtain the numerical results pre- 
sented in SJ5] 

5. Numerical results 

In this section the numerical results on stability and transient growth in Taylor- 
Couette flows are presented. 

5.1. Optimal transient growth in various regimes 
According to the numerical strategy discussed in SJ3]and £4^ the optimal transient growth 



G m ax = sup„ k t G(t) is computed for logarithmically equidistant shear Reynolds numbers 
250 < Re < 2 -'lO 6 and n e {0.2, 0.5, 0.8}. 
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Figure 4. Representative lines in the Rei-Re plane in the case n = 0.5 for which the optimal 
trans ient growth G ma x and corresponding optimal axial wavenumbers fc max are plotted in figures 
5(a) resp. 5(b) 



By studying test cases we find that t £ [0; r ] with tq 



0.85 is 



p J^IT — v and < 

suitable choice to determine the transient growth maximum in time for all considered 
parameter regimes. Optimization in the wavenumbers is carried out by default in the 
range n £ {0, 1 . . . , 8} and k £ [0; 5]. Additionally, as discussed in 1 3.3 the ranges for t, n 



and k are enlarged whenever the optimization terminates near one of the upper bounds. 
By the choice of rotation numbers Rq the linearly stable regimes I, II are parametrized 
considering R n £ [^jlO^] (I) and Rq[-10; -1.1] (II) (see table FT}. Furthermore, 
transient growth is studied in the counter-rotating regime IV nearby the Rei = line 
in the Rei-Re half plane by choosing Rq £ [0.1— -; 0.9^— -]. For a global overview the 

results for R n £ {-3,-1.2,0.8^, 1.2^,3^} are presented. The lines in the Re z - 
i?e -planc defined by this choice for rj = 0.5 are given in figure HI for orientation. Figures 
|5(a)| and |5(b)| show the optimized transient growth G max respectively the corresponding 
optimal axial wavenumber fc max against Re for the considered parameter sets. 

The most prominent feature in the double-logarithmic plots of figure |5(a) are the 
nearly identical asymptotic slopes of the lines in the linearly stable regimes for Rq £ 
{—3,-1.2,1.2,3} showing a characteristic power law G max ~ Re a with a ~ | ± 7% 



(compare dashed line in figure 5(a)). Notably, even in the Ray leigh-unst able counter- 
rotating regime IV (circles in figure 5(a) ) G max seems to approach this scaling as long as 



the computation is not destabilized by dominating linear instability. In fact, for constant 
Re the energy amplifications G max (i?e) in the different regimes differ only by O(l) factors 
and not - as possibly expected - by orders of magnitude. Within the linearly stable regimes 
I and II these deviations are most distinct in the vicinity of the Rayleigh line respectively 
the boundary to regime IV where larger amplifications occur. 

Hence, the numerical results suggest that optimal transient growth in linearly stable 

2 

Taylor-Couette flows roughly follows a common scaling G max ~ Re 3 for Re — > oo. Note 



that this scaling result is in perfect agreement with those by Yecko ( 2004 1 obtained for 
Keplerian flows at fixed i?n 



1.5 in rotating plane Couette geometry. 



5.2. Optimal axial wavenumber 



Beyond the magnitude of transient growth studied in [5.1 the spatial structure of the 
optimal perturbations it max is of great interest. The latter is determined by the optimal 
axial and azimuthal wave numbers fc max , n max which attain the optimal transient growth 
G max shown in figure 5(a) The fc max are plotted in figure |5(b)| with logarithmic hori- 



Transient growth in Taylor-Couette flows 



13 



?7=0.2 



?7=0.2 



77=0.8 





(a) Optimal transient growth G n 



4- 



« 3 



a 2- 



4- 



« 3 



a 2 



oL 



4- 



w 3 



a 2 



oL 













▼ 




i i i? n =-3.0 
T T B n =-1.2 

• • R n =3.2 
J?„=4.8 

* * i? n =12.0 - 


^ 




▼ 










• 



▼ 


""*N 









10* 

Re 
r7 = 0.5 



***** 



A 


A 


-Rn 


= -3.0 


V 


v 


-R<> 


= -1.2 





O 


!.',» 


= 0.8 








R n 


= 1.2 


* 


ft 


R !l 


= 3.0 



— iiih''!' 



******** 

■ imuiiiuiuuim. 



'*»»**«** A A* A AAHmu******** AAA 



Re 
77 = 0.8 




********** 



**"»"hm> a»a>aaaaaaa a 
..>... .>....x.....„,w. 1I „ ii,i ii i 1 111111I i 1 iii 



10 s 



10* 

Re 



io J 



10" 



(b) Optimal axial wavenumber fc n 



Figure 5. Numerical results on (a) optimal transient growth G max and (b) corresponding 
optimal axial wavenumbers fc max against the shear Reynolds number Re for different r\ and 
Rn e {-3,-1.2,0.8^,1.2^,3=^} corresponding to the lines in figure [il in regimes I, If 
and IV; discontinuities in (b) are due to changes in the discrete optimal azimuthal wavenumber 
ra max ; the asymptotic slopes in (a) show a common scaling of G max ~ Re a for a m | for high 
Reynolds numbers Re — > co (dashed line) 



11 



S. Maretzke, B. Hof and M. Avila 







77=0.2 








IS 3D 


-'■■ Solid-body line 
— Rayleigh line 












b 

C 20 

■o 

c 










en 




c 
c 






D 




4 


6 8 1 



Outer Reynolds number Rt t 




Outer Reynolds number R> 



V=0.B 






-■-■ Solid-body line 
--- Rayleigh line 














6 8 1 



Outer Reynolds number Rc^ 



Figure 6. Domain of the quasi-Keplerian regime (II) where the optimal perturbation is axially 
independent (blue shading), i.e. fc max = 0, for r\ — 0.2,0.5,0.8; the boundary line (blue solid 
line) has been determined by a bisecting algorithm with relative accuracy e = 10~ 2 ; in the white 
regions between Rayleigh line (red) and solid body line (black) fc max 7^ 



zontal axes. Note the curves' discontinuities whenever the (discrete) optimal azimuthal 
wavenumbers n max changes. 

The plots reveal a characteristic two-dimensionalization of the optimal perturbations 
in regime II (also observed by Yecko ( 2004[ )): for Re > Re$ the optimal transient growth 
G ma , x (Re) is consistently attained by axially independent perturbations, i.e. fc max = 
(compare Rn = —3 and Rq = —1.2 in figure [5(b) [ ). The transition to fc max = typically 
occurs already for Reynolds numbers as small as Re^ = O(10 3 ). Only near the Rayleigh 
line - that is for —1.2 ^ Rq < — 1 - a sharp divergence of Reo for Rfi —> —1 is observed. 
Here fc max « 1 holds up to the greatest studied shear Reynolds numbers Re = O(10 6 ). 

While fcmax — is only obtained in the quasi-Keplerian regime (II) in regime I (repre- 
sented by Rn = 3, 1.2) fc max seems to (slowly) decay to zero for Re — > 00. At least weak 
axial dependence fc max < 1 is observed for Re > O(10 4 ) in these flows. Once again, the 
asymptotic decay fc max — » is most distinct near the solid body line i?n — > °o and is 
lost near the transition to counter-rotation at Rei — 0. Here an almost constant optimal 
wavenumber fc max = O(l) is observed. 

For further illustration figures [6] and [7] show contour plots of fc max in regimes II resp. I. 
The boundary lines have been computed by a bisecting algorithm with relative accuracy 
e = 10~ 2 . The extent of the shaded regions in figure [6] emphasizes the dominance of 
axially independent perturbations for quasi-Keplerian flows. 

In the counter-rotating regime (IV) we observe a growing fc max with Re. This difference 
might be explained by emerging linear instabilities which first appear for k > 1 in this 
regime and thus render three-dimensional perturbations less dissipative. 

The behaviour of the optimal azimuthal wavenumber n max is not discussed in detail 
here. Notably however, axisymmetric perturbations (corresponding to n = 0) never attain 
significant energy growth G > 0(1) up to high Reynolds numbers Re = O(10 6 ) except for 
a small neighbourhood of Rq = — 1 where the dominant Taylor vortices related instability 
of regime III emerges. On the other hand, for different n^0 usually transient growth 
of the same order is attained. Numerical results indeed suggest that for sufficiently large 
shear Reynolds numbers ?i m a X depends more on the geometrical parameter 77 rather than 
on Re or Rq which parametrize the base flow. In general, an azimuthal wavenumber n 
seems to be optimal if the associated wavelength is A w 27r(n(l — ?7)) _1 ~ 2, i.e. twice 
the gap width, leading to vortices which are of about the same radial and streamwise 
dimension (see e.g. figure [8(a)| centre-right). 

In contrast, the dominant axial wavenumbers k < 1 in regime I correspond to wave- 
lengths of O(10) rather than O(l) gap widths. The optimal perturbations' axial depen- 
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Figure 7. Contour plot of the optimal axial wavenumber fc max attaining optimal transient 
growth G max within the regime I of the Rei-Re parameter space; lines determined by bisection 
at e = 1CP 2 ; discontinuities are due to optimization in the discrete azimuthal wavenumber n 



dence is thus indeed weak compared to azimuthal (and radial) variations. Moreover, the 
stronger the rotational influence on the fluid's stability expressed by smaller 77, the smaller 
are the /c max attained for Re — ¥ 00 (figures 5(b)). The observed two-dimensionalization 



is thus in good agreement with the Taylor-Proudman theorem stating that a rapidly 
rotating inviscid fluid is (preferably) uniform along its rotational axis. 

Changing r\ does not seem to have any further qualitative effects on transient growth 
according to the results in figure^] as long as none of the limits r\ — > {0, 1} is considered. 
A further study of this parameter is therefore omitted in the following. 



5.3. Evolution of optimal perturbations 

In the sequel three different optimal perturbations it max ,i, M m ax.2 and it ma x.3 are con- 
sidered at a constant shear Reynolds number Re = 8000 and 77 = 0.5. The rotation 
numbers are given by Rn.i — —2.0, Rn,2 = 2.0 and Rn,3 = 0.8 corresponding to regimes 
II, I and IV. The optimal wavenumbers are given by fc ma x.i = 0, fe max ,2 ~ 0.464 and 
fcmax,3 ~ 1-200 respectively n maXj i = n max2 — "max, 3 = 3. The time evolution of these 
modes is computed by eigenmode decomposition at a polynomial resolution N = 50. 

In figurelHlthe resulting real parts of it maXj i, M m ax,2 and u maXj 3 are shown at a sequence 
of snapshots tj throughout the transient growth evolution. The flow fields are plotted 
in radial-azimuthal projection (top) and radial-axial projection (bottom) with z on the 
horizontal axis except for w maX) i where the latter is omitted due to the axial indepen- 
dence. The radial-axial plots have been rescaled so that exactly one axial wavelength is 
displayed. Arrow lengths scale with the absolute flow velocities although different scalings 
are applied in figures [8 (a) 1 18 (b) | and |8 (c) The colour map from blue to red marks regions 
with relatively high respectively low energy densities |w m ax,?;| 2 in the current fields. 

The perturbations' total kinetic energy evolutions ||wmax,j(t)|| 2 m relation to the tran- 
sient growth maxima G maXj i are plotted in figure [9| with time scale renormalized by 



2tt 

TO — fle 0.85 (1 



-n) 



The tj considered in figure 



are identified by markers. 



The radial-azimuthal projections in figure 8 reveal essentially similar transient growth 
mechanisms of the considered modes: the optimal initial perturbations have a spiral-like 
structure of 2n streamwise elongated vortices. Recalling the different angular velocities 
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Figure 8. Evolutions of the real part of the optimal perturbations ti ma x,i, it max ,2 and 'Umax, 3 
in radial-azimuthal and radial-axial projection for shear Reynolds number Re = 8000 and radii 
ratio n = 0.5 as examples for transient growth in the regimes II ((a), quasi-Keplerian), I (b) 
and IV ((c), counter-rotating); the subfigures each show subsequent snap shots at times t — tj 
during the transient growth evolution; the tj are also marked in the energy evolution curves 
plotted in figure [9] arrow lengths scale with the flow velocities whereas the colour map from 
blue to red reflects energy densities |u max ,i| 2 ; the radial-axial plots have been rescaled to show 
exactly one axial wavelength in the horizontal axes; computed at N = 50 
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Figure 9. Evolution of the optimal perturbations' kinetic energies ||it ma x,i|| throughout the 
transient growth dynamics for Re = 8000, r\ = 0.5 and Rn,\ = —2.0 (quasi-Keplerian regime 
II), 7?n,2 = 2.0 (regime I) resp. Rn,3 — 0.5 (counter-rotating regime IV); normalization of the 
time axis with to — Re oM$ r (1 _ •, ; snap shot times tj of the velocity fields plotted in figures [8| 
highlighted by markers; computed at N = 50 






Qi, Q of the driving inner and outer cylinders, i.e. O^ > tt > for Rq = —2.0, 
n o > Qi > for Rq = 2.0 respectively Qi > > Ct in the counter-rotating case 
Rn = 0.8 we find that the initial spiral orientations are always misfit to the base flow. 
This "misfit" character is a manifestation of the perturbations' non-modal nature and 
thus typical of transient growth as emphasized by Grossmann ( 2000[ ). The spiral velocity 
fields are tilted by the base flow and thereby gain energy (compare figures M centre-left 
and figure p|. As for the axially- independent perturbation in 8(a) the energy maximum 
then occurs exactly at the turning point of the spirals orientation whereas in cases 2 
and 3 it is attained shortly after this point (centre-right in figure [8]). Subsequently, the 
perturbation is further deformed into a "fit" flow direction, i.e. an eigendirection, and 
meanwhile decays. 

Indeed especially for the two-dimensional perturbation it m ax,i the energy growth and 
decay is a rather sudden phenomenon leading to a sharp peak as depicted in figure 
[9] Possibly this is due to the rapid flow which occurs at the inner cylinder wall (see 
figure 8(a) centre-right) leading to high dissipation around the energy maximum. On the 
other hand, the optimal perturbations w max ,2 and M maX! 3 in the regimes I and IV seem 
to be stabilized in this respect by their axial dependence leading to 40% resp. 115% 
higher growth than the one attained for Rq = —2.0 and slower decay in figure [9] This 
interpretation is stressed by the fact that in spite of the small wavenumber fc max ,2 ~ 0.464 
in case of w ma x,2 up to 86 % of its kinetic energy is temporarily transferred into the axial 
component around the transient growth maximum. 

The characteristic structure of deforming elongated vortices is also reflected in the 
radial-axial projections in figures [8(b)] and 8(c) whereas the elongation occurs spanwise 
in this perspective. A unique feature of the counter-rotating flow (Rn.3 = 0.8) is the 
optimal perturbation's localization near the inner cylinder walls where the base flow 
is locally Rayleigh-unstable. Hence, although the flow remains eigenvalue stable for the 
chosen parameters emerging instabilities already seem to interfere with non-modal growth 
mechanisms. This possibly explains the greater transient growth in regime IV. 
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Figure 10. (a): Numerically computed optimal transient growth G^ ax for axially independent 
perturbations plotted against the shear Reynolds number Re; the curves are independent of Rn 
and parallel for Re ^ O(10 4 ) corresponding to a common scaling G(^ x = a(r/)Re3; (b): Fitted 
scaling coefficients a(r/) for the respective 77 at large Re ^ O(10 4 ); red bars: data determined by 
fits on numerical results for G^^(Re)\ blue line: cubic fit according to ||5.1[) and (|5.2[) 



5.4. Transient Growth Scaling for k = 

The previous numerical results, especially those for the quasi-Keplerian regime II, mo- 
tivate the transient growth a nalys is of axially independent perturbations with k - 
Moreover, it will be shown in [6.2 that G^axj i- e - the optimal transient growth of k 



0. 


disturbance s, is in deed independent of the rotation number Rq . 

In figure 10(a) numerically computed optimal transient growth G^ x is depicted in 
a log-log-plot in the range Re € [250; 8 • 10 6 ] for different 77 e {0.05,0.2,0.5,0.8,0.95}. 
The parallel slopes for high Reynolds numbers Re ^ O(10 4 ) show a common scaling 
G^zZ ~ Re 1 where the proportionality factor may depend but on r\. In order to estimate 



this G^~°(i?e) is computed for logarithmically equidistant Re € [10 5 ;4 • 10 6 

{0.05,0.1,0.15, 

exponents 7(77) 



,0.95}. Fits of the form Gf n =°(i?e) 



for rj € 
a{n){Re) 1 ^ for each 77 yield 



I within errors ^ 0.5 %. Hence, a common exponent 7 = | is faithfully 
assumed and the factor a(rj) is independently determined by another fit. The results are 
plotted in figure |10(b)| where the errorbars give the mean square deviation. 

In order to obtain an analytical formula for G^^.(Re) a third degree polynomial 



(5.1) 



a(rj) = a + a x rj 1 - -r) ) + a 2 r\ 1 



is fitted to the data in figure [10(b) | taking into account the extremum of a at 77 = 1 which 
is due to the system's symmetry with respect to exchanging of r, and r . The result is 



9.218- 10" 



a wy.ZLS- 1U ", ax « 0.1198, and a 2 « -9.072 • 10 -2 (5.2) 

and the corresponding curve is also shown in figure |10(b)| Good agreement between fit 
and data is found especially for 77 > 0.5, possibly due to the lesser impact of the azimuthal 
wavenumber's discreteness on the attainable optimal transient growth compared to 77 < 
0.5. For arbitrary 77 test cases give less than 7% error if the analytical formula is applied 
for Re e [10 4 ; 8 • 10 6 ] and less than 5 % in the interval [10 5 ; 2 • 10 6 ]. 

The maximum amplification of axially independent perturbations G r ^ x = a(j])Re^ 
defines a lower bound for the total (k ^ 0) transient growth G max (Re) in every flow 
regime. Moreover, the estimate can be expected to hold within a factor of O(l) being 
exact in the blue-shaded regions of the quasi-Keplerian regime II in figure [6] 
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6. Analytical results for axially independent perturbations 

The prominent role played by axially independent perturbations together with their 
geometrical simplicity motivates an analytical study ol their properties, which is pursued 
in this section. We begin by applying the conjugated curl operator 

-i* i 

(Vx) r := e - i ^ +kz \Vx)e iinip+k ^ = [ ifc -V \ (6.1) 




to the linearized Navier-Stokes equation (2.5). This eliminates the pressure gradient 
terms, yielding 



(6.2) 



For axially independent perturbations (k = 0) the azimuthal velocity u v is determined 
from u r via the divergence condition 




= V r • u = T> + u r 



iku. 



=o 



-V + U r , 



(6.3) 



and the evolution equation for u r and u z decouple ( Gebhardt fc Grossmann|1993 ) . Using 
C rr — C vv the resulting equations read 






-VC zz u z 



fd t u z 
-Vd t u z 

y (->f + V+lV + )d t U r ) \(-f C rr +V+C rr fV+ + V+C ipr +C rv V+)u r 

The first and second equation, which are equivalent, determine the evolution of u z : 
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(6.5) 



Using the results V + C vr + C rtp T> + = ^V + - f? and [C rr , rd r ] = 2C rr + 2\nA =: 2£° rr 
obtained in |A.1| the evolution equation for u r becomes 
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Further using ft- (rV + rVj r — n 2 ) 
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(6.7) 



This fourth order PDE is supplemented with the boundary conditions u r (ri) — u r (r ) — 
d r u r (ri) = d r u r (r ) = 0, which correspond to the no-slip boundary conditions at the 
cylinders u r (^) = u r (r ) = u v (r t ) = u v (r ) = 0. 
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6.1. Advection of perturbations by the basic flow and universal stability properties 



A remarkable property of the equations (6.5) and (6.7) is revealed by considering the 
transformation u r :— e mAt u r , u z :— e mAt u z . The derivatives then read d r u* ■ 
z mAt (d t + 'mA)u. Jf so substituting into (6.5) and (6.7) yields 



i At 
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(6.8) 
(6.9) 



where f r := (r'D+rD+ — n 2 ) u r and f r :— (rT> + r'D + — n 2 ) u r . As u z and u r satisfy 
equations (6.5) and (6.7) with A = 0, the A-dependence of the perturbation's evolution 
u is entirely described by the factor e~ m . This factor corresponds to a pure advection of 
the perturbation with the shearless, uniformly rotating part of the basic flow v B and thus 
it is locally and globally energy-conserving (|it r | 2 = e mAt e~ mAt \u r \ 2 — |u r | 2 ). Although 
these conclusions might seem obvious at first glance note that they are not true in the 
general three-dimensional case k =/= 0. 

The minor importance of the parameter A has crucial consequences: without loss of 
generality A = can be assumed when analysing the stability of Taylor — Couette flow 
to axially independent perturbations. The remaining parameter B characterizing the 
base flow v B depends only on the shear Reynolds number Re and not on the rotation 



number Rq (see (2.106)) which parametrizes the flow regime. Hence the linear stability of 



Taylor-Couette flow to axially independent perturbations is independent of Rn and thus 
identical in all regimes. Further, the optimal transient growth Gjj^c f° r k — defines a 
lower bound to the absolute maximum G max which is universal in the sense that it only 
depends on r\ and Re. We note that these results can be expected to apply approximately 
also for weakly axially dependent perturbations in the vicinity of k = 0. 

6.2. Global analysis of the evolution equations 
First, consider the evolution of u z described by equation (6.5): C zz is the sum of a self- 



adjoint negative definite operator and a skew-hermitian one. As these do not commute 
C zz is an example of a non-normal operator which nonetheless does not allow for tran- 



sient growth (see appendix A. 2 for a proof). The evolution equation (6.7) for the radial 



component u r may be split into two independent problems 

r 



dtfr = ( £-rr — 29 r - I f r and {rT> + rT> + — n 2 ) u r = f r 



(6.10) 



where the second is of Sturm-Liouville type (the solution is given in appendix A. 4 ) and the 



first resembles equation (6.5 ). Unfortunately, it turns out to be difficult to incorporate the 



boundary conditions for u r in this approach. On the following the discussion of equation 
(6.7) is therefore confined to the limit of asymptotically large Reynolds numbers Re — > oo 



and is studied by means of scale analysis. 

In order to identify and motivate the scales to be studied quantitatively in the WKB- 
Analysis &|7]we consider the energy evolution of a perturbation u = u r e r + u B e„ 



"ip^ip 



dt|H| 2 = 2Re (u, Cu) = -2Re (u, (u ■ V)v B ) + 2Re (u, A r u) (6.11) 

where pressure and convective terms drop out as in the derivation of the Reynolds-Orr 



equation. Using u v — —T>+u r the non-normal term in (6.11) becomes 
N(u) := -2Re(it, (u ■ V)v B ) = lm(u r ,d r u r ) M 



(6.12) 
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while the self-adjoint, dissipative summand reads 

D(u) := 2Re (u, A r u) = -2( n~ 2 (||£>r:D + u r ||H + (n 2 + l)||2? + u r ||^) + ||Z>u r ||JH 

- 4Rc{u r ,r- 1 V + u r ) M + (n 2 + l)\\r- 1 u r \\l). (6.13) 



Assume that u r varies on a typical length scale of order 0{(nRe)~ a ) with a > 0. In 
the limit Re —¥ oo the highest order r-derivative dominates in each term of (6.11) . As 
B ^ Re and ||m|| 2 = ||wr||i[ + n _2 ||rI? + w I .|| 2 I we obtain from equations (6.121 and (6.13) 



N(u) = 7i- 2 0((ni 1 > e) 1 + a ||u r || 2 I ) = 0((ni2e) 1 - ae )||t*|| : 

D(u)= n- 2 0{{nRe) ia \\u r \\l) = 0({nRe) 2a )\\u\\ 2 



(6.14) 



According to (6.13) D(u) is strictly negative, so by virtue of (6.14) dissipation always 



dominates for a> k- On the other hand, the non-normal term N(u) may be positive so 
that for a ^ | growth rates dt In ||m|| 2 = 0((nRe) 1 ~ a ) are possible. 

The question remains how long such growth may last. Let us consider a Fourier- type 
ansatz u r ~ e lmr with wavenumber m = 0((nRe) a ). Note that locally this is valid 
because in the limit Re — ¥ oo boundary effects are confined to thin layers near the 
cylinder walls. Then N(u) is of optimal order in (6.14) and N(u) > if and only if 
n~ x Bm < by virtue of (6.12). 

The total velocity field is u = e 1 ^ nip+kz ^u ~ e 1 ('"i , + mr ) j so the curves of constant phase 
(characteristics) are (locally) given by (p(r) — tp(ri) — n~ 1 m(r — ri). Starting at the inner 
cylinder the set of these lines form streamwise elongated spiral-structures like the vortices 
in figure [8| To attain growth they have to be orientated according to the sign 



sgn(d r <p) 



-sgn^n 



(n m) = sgn(i?) = — sgn(d r il). 



(6.15) 



Thus, the perturbations' characteristics have to be misfit to the base flo w's angular ve- 
locity profile il B 



r~ l v B 



5.3 



Therefore, 



as observed in the numerical computations of { I 
energy amplification may only occur transiently until the perturbation has been sheared 
into the "fit" orientation by advection. Within the advective time scale T = 0(Re~ 1 ), 
i.e. a cylinder rotation period, this shear uniformly distorts the flow profile between inner 
and outer cylinder by a length of order O(l) . Consequently, as the initial streamwise 

) and m — 0((nRe) a ) the time to for the 



elongation of the characteristics is 0(n m 
perturbation to be tilted into fit direction is 



t, 



0,a 



0{n- l mT) = 0{(nRe) c 



(6.16) 



Viscosity prevents transient growth if a > | . Now assume u is an optimal perturbation for 
a < |. Then we can evolve this mode backwards until times of order 0((nRe)~s ) before 
its energy maximum, introduce the result as a new initial condition and thereby attain 
additional growth. Thus, optimal perturbations must vary on length scales 0(nRe)~s 
and to = 0({nRe)~i) is the natural time scale for transient growth. 

Our numerical computations are in perfect agreement with these scaling results. How- 
ever, we cannot obtain an analytical estimate for the optimal transient growth with this 
section's zeroth-order approach. Therefore, in the next section we introduce the time 



scale to to the evolution equation (6.7) and analyse it by means of a first order WKB- 



approximation following the analysis of Chapman (2002, pp. 47-53) for "oblique modes" 
in channel flows. 
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7. WKB- Analysis of axially independent perturbations 

Following the analysis of the previous section we rescale time as i := 5~ 2 t with 6 := 
(nRe)~ s , and rewrite nB := 5~ 3 Bq, where the factor Bq is independent of n and Re (see 



equation 2.106). Substituting these scalings for t and nB in the evolution equation (6.7) 



and multiplying by S 3 yields 



Sdt (r 2 d 2 + 3rd r - (n 2 - 
S 3 id 2 



1)) 



!/.,. 



-d r 

r 



i-Br 



r 2 di + 3rd r 



1)) 



(7.1) 



where we have set A — w.l.o.g. in accordance with |6.1| Note that the highest order 



spatial derivative in (7.1) is now multiplied by the factor S 3 <C 1, which is small in the 



limit of high Reynolds numbers Re — > oo. 

7.1. WKB-ansatz 
We make a WKB-ansatz with amplitude a and rapidly oscillating phase 5~ l (f) 

u r = aexp (<S _1 </>) , (7.2) 

where both a and <fi depend on t and r. Together with the divergence condition this 
yields u v = —T> + u r = 0(d r u r ) = 0(S^ 1 u r ). Hence the scaling a = 5a, with a = 0(1) 
and 4> — 0(1), is required in order that initial perturbations u = u r e r + u, 
energy norm (||it(0)|| 2 = 1). 



v e v have unit 



We now substitute the WKB-ansatz (7.2) into the evolution equation (7.1). Because 
of a, (f> = O(l) the evolution equation needs to be independently satisfied at each order 
in S. At leading order 0(S^ 1 ) the equation reduces to (see Appendix A. 3) 

iBor 



dp, 



LBo 



^> cb(r,t) = cj> (r) 

i - 

By using this solution to eliminate 8~ x terms in ( |7.1[ ) we obtain 
r 2 dt ((d r cj)) 2 a) - r 2 (9 r 0) 4 a = 5 1 (6r(d r c/)) 3 a + Qr 2 {d r 4>) 2 {d 2 r 4>)a 



(7.3) 

4r 2 (d r c/)) 3 (d r a)) 
3r{d r (j))a) (7.4) 



- 5 1 d- t (2r 2 {d r 4>){d r a) + r 2 (d 2 (j>)a 
which reads at next leading order O(8 ) = 0(1) 

(d r cj))dta — (d r (j)) 3 a — 2(djd r (/))a. 
i(d r (j)), dt = i(dtd r (j))d T = -^d T ( |Chapman|[2002| p. 49) yields 
d T a r 3 _ 2 2 t _^ a (r) 



Defining r 



2Bn 



a(r,r) 



exp 



6B n 



(7.5) 



(7.6) 



According to this solution a becomes singular for r — > 0, which may raise doubts about 
its physical correctness. However, in this limit the underlying separation of orders in the 



WKB-approximation breaks down so that 0(S X ) terms in (7.4) or even in the leading 
order equation have to be considered. These bound the blow-up leading to an overall 



nearly singular amplitude behaviour in the complete linearized dynamics given by (6.7). 



In numerical simulations this manifests itself in increasingly sharp peaks of the optimal 
perturbation's energy for Re — > oo as visualized in fig ure [TT] the larger Re the longer 
the blow-up seems to follow the singular WKB-solution ( 7.6 1 before the energy growth is 
capped near the maximum blow-up time to- Most prominently for Re — 1024000 it is only 
in a neighbourhood (1 ± 0.05)io about the maximum that the singularity is smoothed 
out by additional terms resulting in the sharpest peak in figure [TT] 
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Figure 11. Energy blow-up of numerically computed optimal axially invariant perturbations 
for Rn = —2.0, r\ — 0.5 and different shear Reynolds numbers Re; the time axis is normalized by 
the respective energy max imum to; the increasingly sharp peaks reflect the singular behaviour 
of the WKB-solution (|7.6[) except for O(S) neighbourhoods of the maxima 



7.2. Construction of optimal perturbations 



Assume that the amplitude's growth according to equation ( 7.6 1 is capped as soon as the 



next-order terms become relevant. Then the optimal energy growth is attained if 

(a) The blow-up occurs at a common time t over the whole radial domain r € {ri,r ) 

(6) The 0(5 1 ) terms in (7.4 1 are of the highest attainable order in t 

Condition (a) ensures that no averaging effects of the spatial integral evaluated for 

the computation of ||it(£)|| 2 limit the global energy maximum in time. It is equivalent to 

d r (j>{rj ) = d r Mr) + ^*o = 0, so that <j> = ^t + C and w.l.o.g. = ~^(t- to). 
On the other hand, condition (b) implies that the blow-up is capped as late as possible 

in the evolution in r. Let us consider the O(c) 1 ) terms in equation (7.4 1 



- 5 1 d- t (2r 2 (d r (t>)(d r a) + r 2 {d 2 r <j))a + Zr(d r 4>)a) 

2B Q 



5 1 d T ( 2rd r a + (d r r)a + -to) . (7.7) 



Recalling that a = 0(t~ 2 ) and d r a = 0((d r T)r~ 3 ) as r —> we find that the leading 
order terms in ( |7.7[ ) are 0((5 1 (9,-t)t~ 3 ), whereas the left hand side of (7.4) is of order 
r -1 . Hence, the O^S 1 ) terms become significant as soon a,s t — 0({8d r T)?). Accordingly, 
to attain the most sustained blow-up d r T should be as small as possible for r —> 0, i.e. 



0= lim(-i<3 r r) = lim | d 2 

T^0 T^O 



drfo ) = 9> + -drfa. (7.8) 



Equation (7.8) is also satisfied for 



IB 



°t + c. Hence, this is indeed the optimal initial 



phase giving the optimal perturbation according to WKB-theory 
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Note that the boundary conditions are satisfied if and only if a(ri) — a(r ) — d r a(ri) 



(7.9) 



d r a(r ) = so that (7.9) is indeed an approximate solution to the complete boundary 



value problem for t — to = 0(1) if a is suitably chosen. 
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7.3. Boundedness of the blow-up 
According to (7.8) we then have d r r — O(t) for r -> so that the growth is not capped 



before r = O(S). However, it remains to be shown that no further blow-up occurs beyond 
the domain of the WKB-solution ( |7.9| . For times i— i = O(S) we obtain <9™a = 0(5~ 2 ) 
and^exp(-^_(t-io)) 

the scaling St :— t — i , u r 



0{l) for all n € N so that u ri d™u r = O^ -1 ). Therefore 



is employed in equation (7.1) giving 



d~ t (r 2 d 2 + 3rd r - (n 2 - 1)) 



iS 



(r 2 8 2 + 3rd r - (n 2 - 1)) u r + 0(S 3 ). (7.10) 



Setting f r :— (r 2 d 2 + 3rd r — (n 2 — 1)) u r the leading order solution of (7.10) is given 



by f r (r, i) = f r , (r) exp (- ^ (i -i Q )). The operator r 2 d 2 + 3rd r - (n 2 - 1) is of Sturm- 
Liouville type so that a Green's function G(r, r') exists such that 

u r (r,i) = f° G(r,r')/ r ,o(r')exp (-^{i- i )\ r',\r. 



(7.11) 



G is given in Appendix §A.4| Note that with this ansatz only two boundary conditions 
may be satisfied. However, this affects only a thin boundary layer in the vicinity of the 
cylinder walls where significant growth is inhibited already for 0(i — to) = O(l) due to 
the no-slip condition. Thus, the present focus lies on the inner solution in the first place. 
By (7.11) the components u r and u v ~ (1 + rd r )u r are given by i 2 -kernel integral 
operators applied to f r . Consequently, they are L 2 -continuous in f r so that ||m|| 2 depends 
continuously on t. Hence, there is no further blow-up in the time scale i — i — 0(6). 

7.4. A scaling for optimal transient growth 



According to (7.9) the optimal perturbation's components u r and u v ~ (1 + rd r )u r have 
grown to 0(S~ L ) by the optimal (blow-up) time. This yields the optimal transient growth 



G^ = supG(t)=su P ||u(t)|| 



(a) 



/>0 



(>0 



sup(\u r {t)\ 2 



+ \u v (t)\ 2 ) = 0(S- 2 ). 



(7.12) 



by condition (a). Since the WKB-approximation applies for S — > and S = (nRe) 3 it 
has been shown that the optimal transient growth of axially independent perturbations 
scales like G^° ~ i?e 3 in the limit of high Reynolds numbers Re — > 00. This result is in 



perfect agreement with our numerical computations (see ^5.4) 



Notably the scaling exponent a = | is independent of 77 and of Rq (see ^6.1) and 
equal for all azimuthal wavenumbers. In accordance, our numerical results show that as 
Re — > 00 the optimal azimuthal wavenumber n max becomes constant; the asymptotic 
value is selected only by the geometry (specified by 77). 

7.5. Numerical validation 



In order to validate the WKB-solution ( 7.9 1 we compute the initial phases Im(ln u r (r, 0)) 
of numerically determined optimal perturbations u = u r e r + U v e v , as proposed by 
Chapman (2002 p. 51 f.). By equation (7.9) this should yield 



B () 



lm(lnu r (r,0)) =t a + 0(S). 



(7.13) 



Due to the non- uniqueness of the complex logarithm relation (7.13) needs to be assumed 
satisfied for some r € (ri,r ). We choose ro = \(ri + r ). 



In figure 12 the blow-up time to computed from (7.13) is plotted against the radial 



coordinate r (solid red lines). This WKB-prediction is compared for Reynolds numbers 
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Figure 12. Blow-up times to of numerically determined optimal axially independent perturba- 
tions (k — 0) for Rn = —2.0, r\ = 0.5, n = 3 and different shear Reynolds numbers Re; results 
according to the WKB-prediction (7.131 ("by WKB-phases") are contrasted with the numeri- 
cally observed transient growth maximum ( "effective" ) ; "effective ±5" : order of the expected 
error range due to finite Re in the WKB-approximation (S — (nRe)~ ' 3 ) 



Re G {10 3 ,10 4 ,10 5 ,10 6 }, corresponding to 5 £ {0.069, 0.032, 0.015, 0.007}, to the optimal 
time determined numerically from the full equations (dashed blue line). The expected 
error ranges are denoted by [i + 6; to — 5] (dash-dotted blue lines). Excellent agreement 
between the numerical results and WKB-solution within the predicted error of order 
S and convergence for Re — ¥ oo is found. Significant deviations are confined to a 0(S)- 
neighbourhood of the cylinder walls in which growth is prevented a priori by the boundary 
conditions. Hence, the initial phase's behaviour as a key property of the derived WKB- 
approximation has been numerically verified. 



8. Discussion 

Certain Rayleigh-stable Taylor-Couette flows tend to become turbulent at moderate 
Reynolds numbers Re = 0(1000) ( |Taylor|1936"l|Borrero-Echeverry fc Schatz|2010"l |Burin 



& Czarnocki 2012), whereas in case of the quasi-Keplerian regime II the existence of 
turbulence remains debated ( |Ji et aL|2006 Paoletti fc Lathrop|2011 ). At the same time, 
Rayleigh-unstable but eigenvalue stable counter-rotating Taylor-Couette flows are known 
to undergo subcritical transition (Coles 1965). 

In this work the optimal linear transient growth G max , i.e. the maximum non-normal 
energy amplification of infinitesimal perturbations, has been studied as a function of the 
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shear Reynolds number Re, the cylinder radii ratio r\ and the rotation number Rq cover- 
ing all parameter regimes of Taylor-Couette flow. We find that accurate transient growth 
computations are numerically feasible up to Re = O(10 6 ) - even though the character- 
istic Y-shaped eigenvalue spectrum of the linearized Navier-Stokes operator may not be 
resolved for such Reynolds numbers. In contrast, previous studies stress the spectrum's 



significance for transient growth at least in channel flows (Reddy & Henningson 1993). 
In case of Taylor-Couette flow we cannot confirm this close correlation, but find that 
(fortunately) the transient growth maximum G max may long be converged for resolu- 
tions where the approximated spectrum is still far from its natural shape. We exploit 
this numerical circumstance to examine the optimal transient growth for large Re and 
find an asymptotic scaling G max ~ Re a for Re Js O(10 4 ) with a ~ = for any considered 
n £ {0.2,0.5,0.8} and all linearly stable flows. 

This reveals energy growth of the same order in all regimes and allows for arbitrary 
transient amplifications if Re is sufficiently large. Moreover, the optimal perturbations' 
dynamics discussed in j ]5.3| suggest that the underlying growth mechanisms are essentially 
the same in the studied regimes I, II and IV except for additional amplifying effects of the 
Rayleigh instability in the latter. A distinction is found in the optimal axial wavenumbcr 
fcmax which reflects the axial dimensionality of the optimal perturbations attaining max- 
imum energy amplification. While fc max = 0, corresponding to axially invariant modes, 
holds true practically within the whole quasi-Keplerian regime (II) above Re = 0(1000), 
only weakly three-dimensional optimal perturbations < fc max < 1 for Re -> oo are 
found in the likewise Rayleigh-stable regime I. For counter-rotating flows greater fc max 
turn out to attain even higher energy maxima. We conjecture this is due to the emerging 
linear instabilities in this regime which arise for significantly axially dependent perturba- 
tions. On the other hand, it remains unclear why axial dimensionality enhances transient 
growth only in regime I and not for equally Rayleigh-stable quasi-Keplerian flows. 

Overall our numerical results reveal an important role of axially independent perturba- 
tions for transient growth in linearly stable Taylor-Couette flow and so the corresponding 
linearized Navier-Stokes equations have been studied analytically in £J6] and [J7] Firstly, 
the analysis has revealed that transient growth and linear stability are indeed indepen- 
dent of Rq in case k = 0. Then we have shown optimal perturbations to blow up and 
decay within the time scale t = 0((nRe)~%). By introducing this scale to the linearized 
evolution equations an optimal transient growth scaling G^^.(Re) = a(rj)Re^ for axially 
independent perturbations has been derived in the limit Re — \ oo following the channel 



flow WKB-analysis of Chapman ( 2002 ) . The results universally apply for all Rq and thus 



in all flow regimes. For the coefficient a{rf) the semi-empirical formula (5.1 1, (5.2 1 has 
been obtained by a cubic fit on the numerical data. 

The expression G{^ x (i?e) = a(rj)Re^ then provides a universal lower estimate for the 
optimal transient growth of general three-dimensional perturbations. This bound attains 
the optimum in most of regime II according to the numerical results. However, while 
quasi-Keplerian flows thus indeed have the smallest possible energy amplification poten- 
tial, it is nevertheless of the same order as in the other regimes. Temporary amplifications 
of disturbances may promote nonlinear instability if growing modes are consistently fed 
by nonlinear energy redistribution. Hence, by our scaling results such a transient growth 
mediated instability is as likely to exist in quasi-Keplerian flows as in any other regime. 
However, axially independent perturbations are possibly not equally fit to feed nonlin- 
ear instabilities as three-dimensional ones dominating in other regimes, e.g. due to their 
sharper growth and decay. These questions should be addressed in the future. 



Meseguer ( 2002 ) indeed found a strong correlation between the experimentally ob- 
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served nonlinear stability boundary and optimal transient growth G max in counter- 
rotating flows. Following these ideas we estimate the threshold shear Reynolds number 
Rex for subcritical transition in quasi-Keplerian flows using our universal scaling result. 
To this end G max is computed numerically at the subcritical stability boundary of Taylor- 



Couette flow (results not shown) according to measurements by |Mallock (f896), Taylor 

(p36l),|Coiesl " 

Avila & Hofl (|2013|). Not surprisingly the correlation is not as strong as observed by 



Mcseguer (2002) (who considered merely data by Coles (f965)) since transition to turbu- 
lence is no sharp or unambiguously defined phenomenon. Additionally, [Burin fc Czarnocki] 
( 2012[ ) have found their experimental results to depend significantly on the applied end- 
cap configurations where the sensitivity is stronger for wider gaps. Our results indeed 
range from G max ~ 54 to G max ~ 155. If we translate this to shear Reynolds numbers the 
uncertainty roughly agrees with the observed endcap effects. Calculating the mean value 
of all computed threshold amplifications yields an a priori estimate for the threshold 
transient growth in an unknown Taylor-Couette flow setting of G max ,T = 92 ± 26. 

Applying the estimate formula for G max we obtain a threshold Reynolds number of 
Re T = a(rj)-i (880 ± 370) giving for instance Re T = 67000 ± 29000 if r\ = 0.7. In 
case of quasi-Keplerian flows recent expe riments have p roce eded far beyond up to Re = 
O(10 6 ) yielding contradictory results (see Ji et al. ( 2006[ ) and Paoletti & Lathrop (2011 )). 
However, lAvila (2012) has shown the state-of-art Taylor-Couette apparatuses unsuited 



for such measurements due to axial endwall effects. On the other hand, the estimated Rex 
lies still within the range of direct numerical simulations. Hence, such may resolve the 
mystery of turbulence in the quasi-Keplerian regime and perhaps validate the significance 
of G max as a measure for subcritical instability of Taylor-Couette flow in general. 



Appendix A. 

A.l. Calculation of the simplified linearized equations 

In this appendix a few supplementary computations for the derivation of the evolution 
equations in section [6] are presented. 

Firstly, the commutator relation [rd r , C rr ] = £° r is shown. Setting a :— n 2 — I + inB 
we obtain 



[l^rr ; ^^VJ 



\v v n2 ~ l P 


in B a 

v v ,rd r 

r v 

= [d 2 r ,rd r ] 




\d r + -) y,rd r 

\ r J r z 


+ 


-d r ,rd r 
r 


- 


~^,rd r 



2dl + 2d r - 2^r = 2C„ + 2mA = 2£° 



(Al) 



Moreover, the expression T> + C vr + C rip T> + can be simplified by 

' 2in 



LJ+C^r + L, ri pD+ — i?4 



-2A 



, 2B 2in\ 2B 4m 

2A + — \V+= —V + - 



(A 2) 
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Lastly, the equality d r - (rV + rV + - n 2 ) = 2C° rr rV^ 



im<M-D+ 



d r - (rV+rV+ - n 2 ) = d r - (r 2 d 2 + 3rd r - (r 



2£° rr rV + + irn -^V 



2B 



2rdl + 28 2 + 6d 2 



1)) 



2(n 2 - 1) 



= 2rd 3 r +8d 2 - 2[n2 1} U - 



Ain\ 



I = 2 ( dt + -d r 

r 



n 



r 

2 -l 



holds since 



2(n 2 



inB 



2inB 



4n 2 



rV 4 



= 2 &, 



1 n 2 - 1 

-O r j 



r 



= 2rff + 8df - 



2 2(n 2 -l)A„ I 



{rd r + 1) 
9 r — 



An 2 



(A 3a) 



(A3b) 



A. 2. Analysis of the axial evolution equation 

Consider the operator C zz and the axial component u z from the evolution equation 
(6.5) on the Hilbert space H introduced in [2.1 and let u,v £ M0^ 2 {(ri\r )) satisfy 

2 

homogeneous Dirichlct boundary conditions. Define A\ := 'D+V, A2 '■= — \ and B := 
— — v^. A2 and B multiply by a real and strictly negative resp. purely imaginary function. 
Hence, Ai is self-adjoint negative-definite and B is skew-hermitian. For A\ we have by 
partial integration 



(u,Aiv) 



ru*(d 2 + r 1 d r )vdr 



P .i. 



(d r u*)(d r v)rdr 



p.i. 



(rd r u* + d r u*)vdr = (AiU, v)% 



(A4o) 

(A46) 



Equation (A 46 1 reveals Ai to be self-adjoint and for u — v (A 4a I shows its negative- 
definiteness. Thus, C zz is the sum of a self-adjoint strictly negative operator A := A1+A2 
and a skew-hermitian one, B. For the commutator [•, •] we have 



[AB] = {d 2 r 



l d r ) 



jtO. 



Accordingly the adjoint operator C* zz satisfies 

[£*„,£„] = [A - B,A + B] = 2 [A,B] ^ 



(A 5) 



(A 6) 



so that C zz is a non-normal operator. By definition of ||u z ||jfj[ is equal to radial components 
contribution to the total kinetic energy of u. Due to the evolution d t u z — C zz u z we have 

<9 t ||-u 2 ||H = 2Re (u z ,£zzu z ) m = 2Re (u z ,Au z ) u +2Re (u z , Bu z ) u < (A 7) 

v v ' v v ' 

<o =0 

where 2Re(cc,Ta;) = (x,Tx) + (Tx,x) — (x,Tx) — (x,Tx) — for T - skew-hermitian 
has been used. By relation ( A 7 ) there is no transient growth but only monotonic decay 
in the axial component of k = perturbations as claimed in §6.2| 
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A. 3. WKB-equations for the radial evolution equation 
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In the sequel we derive of the WKB-equations (|7.3| and (|7.4|) 

„2fl2 



Application of the operator lrV + r'D + — n 



d 2 + 3rd r - {n 2 - 1)) to the WKB- 



ansatz u r = Saexp ((5 l (f) of S 7.1 yields 

exp(-(T V) (r 2 d 2 + 3rd r - (n 2 - 1)) u r 

= 5- 1 r 2 (a r ^) 2 a + (5° (2r 2 (a r 0)(9 r a) + r 2 (d 2 <f>)a + 3r(d r cj))a) 

+ S 1 (r 2 d 2 + 3rd r - (n 2 - 1)) a + 0{8 2 ) 



The left hand side of equation (7.1 1 thus reads: 



5exp(-5- l (f>)dt (r 2 d 2 + 3rd r - (n 2 - 1)) u r 

= 6- 1 {di<j>) exp(-(5- 1 0) (r 2 d 2 + ird r - {n 2 - 1)) u r + 8 a r 2 d- t ({d r <j>) 2 A) + 

5 x d- t (2r 2 (d r <P)(d r A) + r 2 (d 2 <p)A + +3r(d r <j>)A) + 0{8 2 ) 

and the right hand side: 



exp(-<TV) (S 3 ld 2 
= exp(— J -1 



1 n 2 — 1 
di d r ^ — 



iBo 



} (r 2 2 + ird r - in 2 - 1)) 



iB 



(r 2 d 2 + 3rd r - (n 2 - 1)) u r + S°r 2 (d r 



{dtMr^drtYa 



+ 5 1 (6r(d r cf>) 3 a + 6r 2 (<9 r 0) 2 (<9 2 0)a + 4r 2 (<9 r 0) 3 (d r a)) + 0(5 2 ) 

(A I 



Hence, for the leading order terms = 0{8 ) equation (7.3) is obtained 



iB Q 



r z {d r (j)Ya) 



Or 



i_Bo 



The next order = 0(5°) equation reads 

r 2 0t {{d r 4>) 2 a) - r 2 (d r ^a = - ( drf + ^ ) cxp(-<T V) (r 2 d 2 + 3rd r - (n 2 - 1)) u r 

+ S 1 (6r(d r <p) 3 a + 6r 2 (d r <p) 2 (d 2 (j))a + Ar 2 (d r <p) 3 (d r a)) 

- 5 1 d i (2r 2 (d r c/))(d r a) + r 2 {d 2 r 4>)a + 3r(d r <j>)a) (A 9) 



Consequently, applying di4> + ^r = from expression (|A8[) to equation ( A 9 1 the next 



to-leading order WKB-equation (7.4) follows 



A. 4. Green's function for the radial evolution equation 
In this appendix the Green's function G used in §7.1| is derived and thereby the regularity 
of the approximative solution u r defined by equation (7.11 ) is proven. 

Consider the eigenvalue problem (r 2 d 2 + 3rd r — (ri 2 — 1)) i/j\(r) — —\ip\(r) in the 
interval r G (ri;r ). With p := r 3 , q := (n 2 — l)r, w := r and boundary conditions 
^xifi) = V^aC^o) = this is a Sturm-Liouville problem of the form 



- d r {p ■ (d r tpx)) +q = Xwipx 



(A 10) 



The eigenvalues {A m } meN are thus discrete and the corresponding normalized eigen- 
functions form a complete orthonormal set {ip m } m€ fi with respect to the inner product 
(iph V , iti)h ~ I " Tpf^mwdr of the Hilbert space H introduced in { 2.1 
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A solution to the inhomogeneous problem (r 2 d 2 + 3rd r — (n 2 — 1)) ip = g, ip(ri) = 
ip(r ) = is consequently given by 

vb(r) = f° G(r,r') 5 (r')r'dr' with G{r,r') := - ^ ^n{r')*^ m (r) (An) 

where G is the Green's function. By definition G is continuous and thus bounded on 
[ri;r ] 2 . For the given problem the normalized solution to the eigenvalue problem reads 



7r 2 m 2 \f— 2Tnrj . ( irm 



r 



X m = n 2 - " — and ip m (r) = - — :l ^ 1 ^- sin ( -f^ In — I , meN (A 12) 

in 77 r \ hit] ri 
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